
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>

#include "fftw3.h"
#include "dpi_algo.h"
#include "pub_lib.h"
#include "image.h"

extern cpu_image effect_data;

fft_data_struc* init_fft_data(int w, int h)
{
	float* image_data; //image_data[h][w]
	fftwf_complex* out;
	void* plan = fft_init(w, h, &image_data, &out);

	float* spectrum = malloc(w * h * sizeof(float));
	BYTE* s_norm = malloc(w * h);

	BYTE* image = malloc(w * h);
	fft_data_struc* ret = malloc(sizeof(fft_data_struc));
	ret->w = w;
	ret->h = h;
	ret->plan = plan;
	ret->image_data = image_data;
	ret->fft_out = out;

	ret->image = image;
	ret->spectrum = spectrum;
	ret->s_norm = s_norm;

	return ret;
}

void destroy_fft_data(fft_data_struc* ptr)
{
	fft_exit(ptr->plan);

	free(ptr->image);
	free(ptr->spectrum);
	free(ptr->s_norm);

	free(ptr);
}

//internal calc and draw spectrum
static void id_spectrum(fft_data_struc* p, int pos, int ts, int nw, int nh)
{
	int w, h;
	w = p->w;
	h = p->h;

	fft_r2c(p->plan, p->image_data, p->fft_out);
	calc_fft_spectrum(p->fft_out, p->spectrum, w, h);
	log_float_data(p->spectrum, p->s_norm, w * h);

	if(ts) {
		//binary image
		int i;
		for(i = 0; i < w * h; i++) {
			if(p->s_norm[i] < ts)
				p->s_norm[i] = 0;
		}
	}

	if(w == nw) {
		draw_right_data(p->s_norm, pos);
	} else {
		show_scaled_data(p->s_norm, pos, nw, nh);
	}
}

/*
 * in:
 *	dib: image handle
 *	center:
 *		1: move spectrum center to [w/2, h/2]
 *		0: keep original
 *	pos: right window pos from 1 to 4
 *  ts: threshold
 *		if ts = 0, ignore it
 *	nw, nh: show as nw * nh
 *
 *	out:
 *	 use destroy_fft_data() to free it
 *	 * */
fft_data_struc* draw_spectrum(FIBITMAP* dib, int center, int pos, int ts,
	int nw, int nh)
{
	int w, h;
	w = FreeImage_GetWidth(dib);
	h = FreeImage_GetHeight(dib);

	fft_data_struc* p = init_fft_data(w, h);
	if(center)
		init_image_center_data(dib, p->image_data);
	else
		init_image_data(dib, p->image_data);
	id_spectrum(p, pos, ts, nw, nh);
	return p;
}

/*
 * in:
 *	data at p->image_data
 *	center:
 *		1: move spectrum center to [w/2, h/2]
 *		0: keep original
 *	pos: right window pos from 1 to 4
 *  ts: threshold
 *		if ts = 0, ignore it
 *	nw, nh: show as nw * nh
 *
 *	out:
 *	 use destroy_fft_data() to free it
 *	 * */
void calc_draw_spectrum(fft_data_struc* p, int center, int pos, int ts, 
	int nw, int nh)
{
	if(center)
		move_data_to_center(p->image_data, p->image_data, p->w, p->h);
	id_spectrum(p, pos, ts, nw, nh);
}

/* filter image with BandReject filter
 * r: cutoff frequency
 * wb: width of the band
 * order: Butterworth order
 *
 * out:
 *	pfft->image_data: original data
 *	pfft->image: scaled data
 */
void br_filter_image(fft_data_struc* pfft, FIBITMAP* dib,
	freq_filter_type type, float r, float wb, int order, int reverse)
{
	int w, h, p, q;
	w = pfft->w;
	h = pfft->h;
	filter_freq_newsize(w, h, w, h, &p, &q);

	float scale = (float)p / w;
	r = r * scale;
	wb = wb * scale;

	fftwf_complex* filter;
	if(reverse)
		filter = bpfilter(r, wb, p, q, type, order);
	else
		filter = brfilter(r, wb, p, q, type, order);

	//image filter
	init_image_datab(dib, pfft->image);
	filter_freq_domain(pfft->image, w, h, filter, p, q, pfft->image_data);
	fftwf_free(filter);

	//scaled
	scale_float_data(pfft->image_data, pfft->image, w * h);
}

/*
 * filter image of dib, result at pfft->image
 * this funciton call fftwf_free(filter) to release filter
 * */
void common_filter_image(fft_data_struc* pfft, FIBITMAP* dib, fftwf_complex* filter)
{
	int w, h;
	w = pfft->w;
	h = pfft->h;

	//image filter
	init_image_datab(dib, pfft->image);
	filter_freq_domain(pfft->image, w, h, filter, w, h, pfft->image_data);
	fftwf_free(filter);

	//scaled
	scale_float_data(pfft->image_data, pfft->image, w * h);
}

